# # Cardiac mechanics Benchmark - Problem 3
# This code implements problem 3 in the Cardiac Mechanic Benchmark paper
#
# > Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S,
# Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software:
# benchmark problems and solutions for testing active and passive material
# behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.
#
import dolfin
try:
from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
from dolfin import DirichletBC, Constant, interpolate, Mesh
import pulse
from fenics_plotly import plot
geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])
# geometry = pulse.geometries.benchmark_ellipsoid_geometry()
2021-11-23 08:26:24,059 [625] INFO pulse.geometry_utils:
Load mesh from h5
# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 2.0
material_parameters["bf"] = 8.0
material_parameters["bfs"] = 4.0
material_parameters["bt"] = 2.0
activation = Constant(0.0)
material = pulse.Guccione(
parameters=material_parameters,
active_model=pulse.ActiveModels.active_stress,
activation=activation,
)
# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
return DirichletBC(
V,
Constant((0.0, 0.0, 0.0)),
geometry.ffun,
geometry.markers["BASE"][0],
)
# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve problem
pulse.iterate.iterate(problem, lvp, 15.0, initial_number_of_steps=50)
pulse.iterate.iterate(problem, activation, 60.0, initial_number_of_steps=50)
2021-11-23 08:26:24,330 [625] INFO pulse.iterate: Iterating to target control (f_28)...
2021-11-23 08:26:24,333 [625] INFO pulse.iterate: Current control: f_28 = 0.000
2021-11-23 08:26:24,334 [625] INFO pulse.iterate: Target: 15.000
2021-11-23 08:27:14,203 [625] INFO pulse.iterate: Iterating to target control (f_20)...
2021-11-23 08:27:14,204 [625] INFO pulse.iterate: Current control: f_20 = 0.000
2021-11-23 08:27:14,205 [625] INFO pulse.iterate: Target: 60.000
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 518),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 553),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 580),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 615),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 650),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 685),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 728),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 779),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 853),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 904)],
[Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 516),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 551),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 578),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 613),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 648),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 683),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 726),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 777),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 851),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 902)])
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
endo_apex_marker = geometry.markers["APEX_ENDO"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)
print(
("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
endo_apex_pos[0],
),
)
Get longitudinal position of endocardial apex: 19.3891 mm
epi_apex_marker = geometry.markers["APEX_EPI"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)
print(
("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
epi_apex_pos[0],
),
)
Get longitudinal position of epicardial apex: 20.7619 mm
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()